Cores

This is figure 3 in the main text

# make Figure 3 Statistical results of cores

library(dplyr)
library(readr)
library(knitr)
library(ggplot2)
library(reshape2)
library(readxl)
library(fitdistrplus)
library(tidyr)
library(classInt)
library(cowplot)
library(gridExtra)

# read in the data

file_name <-
  here::here("analysis/data/raw_data/artefacts-after Sam_Dec2018_quina_modified-2019.xls")

flakes <- as.data.frame(read_excel(file_name, sheet = "flake basics"))
cores <- as.data.frame(read_excel(file_name, sheet = "core basics"))
debris <- as.data.frame(read_excel(file_name, sheet = "chunk&debris"))
retouch <- as.data.frame(read_excel(file_name, sheet = "retouch"))
core_flake_scars_1 <- as.data.frame(read_excel(file_name, sheet = "core flake scars"))
edge_angle_1 <- as.data.frame(read_excel(file_name, sheet = "edge angle"))
retouched_pebble_chunk <- as.data.frame(read_excel(file_name, sheet = "Pebble&chunk"))

scar_dir <- as.data.frame(read_excel(file_name, sheet = 'scar direction'))

core_flake_scars <- core_flake_scars_1[core_flake_scars_1$number %in% cores$number, ]

edge_angle <- edge_angle_1[edge_angle_1$number %in% retouch$number, ]


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# basic details of the dataset

artefact_ids <- list(flake_ids = flakes$number,
                     core_ids = cores$number,
                     debris_ids = debris$number)

retouch_ids <- retouch$number

# how many unique artefacts?

totals_of_each_type <- lapply(artefact_ids, function(i) length(i))

count_unique_artefacts <- sum(unlist(totals_of_each_type))

# debris pieces that are retouch pieces

debris_with_retouch <- intersect(retouch_ids, artefact_ids$debris_ids)

proportion_debris_retouched <- length(debris_with_retouch) / count_unique_artefacts

# proportion of each type

type_table <- data.frame(type = names(as.data.frame(totals_of_each_type)),
                         count = as.numeric(as.data.frame(totals_of_each_type)),
                         proportion = round(as.numeric(prop.table(as.data.frame(totals_of_each_type))),3))

# add total row at the bottom

type_table$type <- as.character(type_table$type)

type_table <- rbind(type_table, c("total", count_unique_artefacts, round(sum(type_table$proportion),2)))

# print pretty table

# print("table of lithic types")

# kable(type_table)

### raw materials

raw_material_table <- data.frame("category" = character(),
                                 'basalt' = numeric(),
                                 'chert' = numeric(),
                                 'limestone' = numeric(),
                                 'quartz' = numeric(),
                                 'sandstone' = numeric(),
                                 stringsAsFactors = FALSE)


res <- flakes %>%
  mutate(complet_flk = ifelse(type1 %in% c("flake", "blade", "KBW", "tablet", "leva", "leva?",
                                           "DBD", "crest"), "complete flake",
                              ifelse(type1 %in% c("flake brk", "blade brk", "LVFB", "proximal"),
                                     "flake breaks", "other"))) %>%
  dplyr :: count(material, complet_flk) %>%
  arrange(complet_flk) %>%
  group_by(complet_flk) %>%
  spread(material, n) %>%
  replace(., is.na(.), 0) %>%
  rename(category = complet_flk)

raw_material_table <- bind_rows(raw_material_table, res[1:2,])

res <- cores %>% filter(!is.na(material)) %>%
  dplyr :: count(material) %>%
  spread(material, n) %>%
  mutate(category = "cores")

raw_material_table <- bind_rows(raw_material_table, res)

res <- debris %>% filter(!is.na(material)) %>%
  dplyr :: count(material) %>%
  spread(material, n) %>%
  mutate(category = "debris")

raw_material_table <- bind_rows(raw_material_table, res)

res <- flakes %>%
  filter(type1 %in% c("ret", "ret brk", "BUR", "LVT", "ret leva?", "ret pro-lepoin")) %>%
  dplyr:: select(material, number) %>%
  dplyr :: count(material) %>%
  spread(material, n) %>%
  mutate(category = "retouched flakes and breaks")

raw_material_table <- bind_rows(raw_material_table, res)

res <- retouched_pebble_chunk %>%
dplyr :: count(material) %>%
  spread(material, n) %>%
  mutate(category = "retouched chunks")

raw_material_table <- bind_rows(raw_material_table, res)

###tools raw materials

res <- flakes %>%
  mutate(type2 = replace(type2, is.na(type2),"unidentifiable" )) %>%
  filter(type1 %in% c("ret", "ret brk", "BUR", "LVT", "ret leva?", "ret pro-lepoin")) %>%
  mutate(tooltype = ifelse(type2 == "backed knife", "backed knife",
                      ifelse(type2 %in% c("bec", "bec,scp"), "bec",
                       ifelse(type2 %in% c("borer","borer,dent", "borer,dent,scp",
                                           "borer,end,scp","borer,notch",
                                           "borer,notch,scp","borer,scp",
                                           "borer,scp,end","borer,strt" ),
                              "borer",
                        ifelse(type2 %in% c("BUR", "BUR? CORE?","BURIN", "burin"),
                               "bur",
                         ifelse(type2 %in% c("cp", "chopper"), "chopper",
                         ifelse(type2 =="cleaver", "cleaver",
                          ifelse(type2 %in% c("dent", "dent,borer,scp","dent,end","dent,notch",  "dent,scp",
                                              "dent/NOT", "dent/POINT"  ),
                                 "denticulate",
                                 ifelse(type2 %in% c("end", "end,borer" ,
                                                     "end,borer,scp" ,"end,scp" ,
                                                     "endscp" ), "endscraper",
                                ifelse(type2 %in% c("leva point","point" ,"POINT"
                                                    , "point break",
                                                    "point,broken", "point?"  ),
                                       "point",
                                       ifelse(type2 %in% c("moon","scp",
                                                           "sidescp" , "sidescp-cvx"  ),
                                              "scraper",
                                   ifelse(type2 %in% c("NATBACK", "nature backed"),
                                          "natural backed",
                                          ifelse(type2 %in% c( "notcch", "notch" , "notch,borer"  , "notch,end","notch,scp" ,
                                             "notch/DENT"),
                                             "notch",
                                             ifelse(type2 %in% c("TANG", "TANGED",  "TANGED POINT", "TANGED POINT /NOT",
                                               "tanged point?", "TANGED POINT?"),
                                               "tanged point", type2)))))))))))))) %>%
  dplyr :: count(material, tooltype) %>%
  spread(material, n) %>%
  replace(., is.na(.), 0) %>%
  rename(category = tooltype)

raw_material_table <- bind_rows(raw_material_table, res)

raw_material_table[is.na(raw_material_table)] <- 0 ##convert na to 0

raw_material_table$total <- rowSums(raw_material_table[,-1])

x <- raw_material_table

prop_raw_material <- cbind(category = x[, 1], round(x[, -1]/rowSums(x[, 2:(ncol(x)-1)])*100, digits = 1))

# for a table like what we have in the paper, SI Table S1, we
# need to do this

# put rows in a specific order
row_order <- c("cores", "complete flake", "flake breaks", "debris", "retouched chunks",
               "retouched flakes and breaks",
               "baked knife", "bec", "borer", "bur", "chopper", "cleaver", "denticulate",
               "endscraper", "natural backed", "notch", "point", "scraper", "tanged point", "unidentifiable")

# for counts
raw_material_table_counts_tbl <-
raw_material_table %>%
  mutate(category = factor(category, levels = row_order)) %>%
  arrange(category) %>%
  filter(!is.na(category)) %>%
  mutate(others = quartz + sandstone) %>%
  dplyr::select(category, chert, basalt, limestone, others, total)

# for proportions
raw_material_table_props_tbl <-
prop_raw_material %>%
  mutate(category = factor(category, levels = row_order)) %>%
  arrange(category) %>%
  filter(!is.na(category)) %>%
  mutate(others = quartz + sandstone) %>%
  dplyr::select(category, chert, basalt, limestone, others, total)

# combine both to produce Table S1
raw_material_table_props_and_counts_tbl <-
tibble(
  type = raw_material_table_counts_tbl$category,
  chert_count = raw_material_table_counts_tbl$chert,
  chert_perc = raw_material_table_props_tbl$chert,
  basalt_count = raw_material_table_counts_tbl$basalt,
  basalt_perc = raw_material_table_props_tbl$basalt,
  limestone_count = raw_material_table_counts_tbl$limestone,
  limestone_perc = raw_material_table_props_tbl$limestone,
  others_count = raw_material_table_counts_tbl$others,
  others_perc = raw_material_table_props_tbl$others,
  total = raw_material_table_counts_tbl$total
)


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### cores attribute summary
core_summary <- cores[0, c(6:17)]

core_summary["mean",] <- c("mean" = apply(na.omit(cores[,c(6:17)]), 2, mean))

core_summary["sd",] <- c("sd" = apply(na.omit(cores[,c(6:17)]), 2, sd))

core_summary["CV",] <- core_summary["sd",]/core_summary["mean",]

core_summary <- rbind(core_summary,
                      sapply(na.omit(cores[,6:17]), function(x) quantile(x, c(0.25, 0.5, 0.75))))


########################################################################



# function to compute frequencies and make a plot
f <-  function(the_data, the_column) {
  data.frame(table(the_data[the_column])) %>%
    ggplot(aes(x = reorder(Var1, Freq), y = Freq)) +
    geom_bar(stat = "identity") +
    ylab("count") +
    xlab(the_column) +
    theme_bw() +
    theme(axis.text.x = element_text(size = 20, angle = 45, hjust = 1, vjust = 1),
          axis.text.y = element_text(size = 20),
          axis.title.x = element_text(size = 25),
          axis.title.y = element_text(size = 25)
    )
}

##################################################################################
# These lines will make figure for cores


# some plots of frequencies for cores

## geometry

p_core_geometry <- f(cores, "core_geometry") + xlab("Core geometry") + ylab("Count")

core_geom_table <-  table(na.omit(cores$core_geometry))

## platform number

p_core_platform <- ggplot(cores, aes(platform)) +
  geom_bar() + theme_bw() +
  xlab("Number of platform") +
  ylab("Count") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))


### cortex proportion of cores

p_core_cortex <- ggplot(cores, aes(cortex_percentage)) +
  geom_histogram(binwidth = 10) + theme_bw() +
  xlab("Cortex proportion (%)") +
  ylab("Number of cores") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))

core_cortex_table <-  table(na.omit(cores$cortex_percentage))

## platform type

cores <-
  cores %>%
  mutate(platform_preparation = case_when(
    platform_preparation == "over" ~ 'overhang\nremoval',
    platform_preparation == "plane" ~ 'plain',
    platform_preparation == "faceted" ~ 'facetted',
    platform_preparation == "cortex" ~ 'cortical',
    platform_preparation == "dehid" ~ 'dihedral',
    TRUE ~ as.character(platform_preparation)
  )) %>%
  filter(platform_preparation != "n") %>%
  drop_na(platform_preparation)

# "over"    "plain"   "plane"   "faceted"  "n"    "cortex"  "dehid".
# 'overhang removal', 'plain' (same as plane?), 'facetted', 'cortical', 'dihedral'

p_core_plat_prep <- 
  data.frame(table(cores["platform_preparation"])) %>%
    ggplot(aes(x = reorder(platform_preparation, Freq), y = Freq)) +
    geom_bar(stat = "identity") +
    theme_bw() +
  xlab("Platform preparation") +
  ylab("Count") +
  coord_flip()

core_plat_type <-  table(na.omit(cores$platform_preparation))

## platform rotation

p_core_rotation <- ggplot(cores, aes(rotations)) +
  geom_bar() + theme_bw() +
  xlab("Number of rotation") +
  ylab("Number of cores") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))

core_plat_rotat <-  table(na.omit(cores$rotations))

## core reduction type

p_core_reduction <- f(cores, "type") + xlab("Reduction type") + ylab("Count")

core_type_table <- table(na.omit(cores$type))

## core scar number

p_core_scar_number <- ggplot(cores, aes(scar_number)) +
  geom_bar() + theme_bw() +
  xlab("Scar number on cores") +
  ylab("Count") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))

core_scar_table <- table(na.omit(cores$scar_number))

## core flake scar length

core_flake_scars_melt <- melt(core_flake_scars_1)

colnames(core_flake_scars_melt)[3] <- "scar_length"

p_core_scar_length <- ggplot(core_flake_scars_melt[core_flake_scars_melt$scar_length < 500, ],
       aes(scar_length)) +
  geom_density() + theme_bw() +
  xlab("Length of scar on cores (mm)") +
  ylab("Density") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))

# Figrur 3 summary of cores is made here
figure_cores <-
  cowplot::plot_grid(p_core_platform,
                     p_core_scar_number,
                     p_core_scar_length,
                     p_core_cortex,
                     p_core_plat_prep,
                     labels = "auto",
                     label_size = 26,
                     ncol = 3,
                     align = 'v',
                     axis = 'lr',
                     #rel_heights = c(1, 1.5),
                     vjust = 1.5,
                     hjust = -0.5,
                     scale = 0.95) #

 # Fig 3 Statistical results of cores

save_plot(here::here("analysis/figures/fig-3-cores.png"), figure_cores,
          base_aspect_ratio = 1,
          base_height = 13,
          base_width = 21)

knitr::include_graphics(here::here("analysis/figures/fig-3-cores.png"))

Flakes

This is figure 9 in the main text (statistical results for flakes).

library(readxl)
library(tidyverse)

file_name <- here::here("analysis/data/raw_data/Copy of artefacts original.xls") 

flakes <- read_excel(file_name, sheet = "flake basics")
cores <- read_excel(file_name, sheet = "core basics")
debris <- read_excel(file_name, sheet = "chunk&debris")
retouch <- read_excel(file_name, sheet = "retouch")
# Figure 9

flakes_clean <- 
  flakes %>% 
  mutate(mass = parse_number(mass),
         `Length/thickness ratio` = length / `Thickness at 50% max dim`)

flake_hist_mass <- 
ggplot(flakes_clean,
       aes(mass)) +
  geom_histogram() +
  theme_minimal() +
  xlab("Mass (g)") +
  ylab("Number of flakes")

flake_dens_max <- 
ggplot(flakes_clean,
       aes(`max dimension`)) +
  geom_density() +
  theme_minimal() +
  xlab("Maximum dimension (mm)") +
  ylab("Density")

flake_lw_ratio_median <- 
  median(flakes_clean$`Length/thickness ratio`,
         na.rm = TRUE)

flake_lw_ratio <- 
ggplot(flakes_clean,
       aes(`Length/thickness ratio`)) +
  geom_histogram() +
  theme_minimal() +
  geom_vline(xintercept = flake_lw_ratio_median, 
             colour = "red") +
  xlab("Length/thickness ratio") +
  ylab("Number of flakes")

flakes_clean_long_thick <- 
  flakes_clean %>% 
  dplyr::select(starts_with("Thickness")) %>% 
  pivot_longer(everything()) %>% 
  mutate(name = str_remove(name, "Thickness at ")) %>% 
  mutate(name = str_remove(name, " max dim"))

flake_thick_dens <- 
ggplot(flakes_clean_long_thick,
       aes(value,
           colour = name)) +
  geom_density() +
  theme_minimal() +
  xlab("Thickness (mm)") +
  theme(legend.position = c(0.75, 0.5)) +
  ylab("Density") +
   guides(colour = guide_legend("Thickness at % of\nmaximum\ndimension"))

flakes_clean_long_width <- 
  flakes_clean %>% 
   dplyr::select(starts_with("Width")) %>% 
  pivot_longer(everything()) %>% 
  mutate(name = str_remove(name, "Width at ")) %>% 
  mutate(name = str_remove(name, " max dim"))

flake_width_dens <- 
ggplot(flakes_clean_long_width,
       aes(value,
           colour = name)) +
  geom_density() +
  theme_minimal() +
  xlab("Width (mm)") +
  scale_x_continuous(limits = c(0, 150)) +
  theme(legend.position = c(0.75, 0.5)) +
  guides(colour = guide_legend("Width at % of\nmaximum\ndimension")) +
  ylab("Density")

flake_cortex_prop <- 
  ggplot(flakes_clean,
       aes(`cortex percentage`)) +
  geom_bar() +
  theme_minimal() +
  xlab("Cortex proportion (%)") +
  ylab("Number of flakes")

flakes_clean_facetted_plat <- 
  flakes_clean %>% 
  mutate(platform = case_when(
    str_detect(platform, "fac") ~ "facetted",
    .default = "unfacetted"
  ))

flake_dens_max_fac <- 
ggplot(flakes_clean_facetted_plat,
       aes(`max dimension`,
           colour = platform)) +
  geom_density() +
  theme_minimal() +
  xlab("Maximum dimension (mm)")  +
  theme(legend.position = c(0.75, 0.5)) +
  ylab("Density") 

flake_scar_count <- 
  ggplot(flakes_clean,
       aes(`scar number`)) +
  geom_bar() +
  theme_minimal() +
  xlab("Number of scars") +
  ylab("Number of flakes")

flakes_clean_5_scars <- 
flakes_clean %>% 
  mutate(more_than_five_scars = ifelse(`scar number` <5, 
                                       "less than 5",
                                       "more than 5")) %>% 
  drop_na(more_than_five_scars)

integer_breaks <- function(n = 5, ...) {
  fxn <- function(x) {
    breaks <- floor(pretty(x, n, ...))
    names(breaks) <- attr(breaks, "labels")
    breaks
  }
  return(fxn)
}

flake_scar_by_mass <- 
  ggplot(flakes_clean_5_scars,
       aes(`scar number`, 
           group = `scar number`,
           mass
           )) +
  geom_boxplot() +
  theme_minimal() +
  xlab("Number of scars") +
  ylab("Mass (g)") +
  scale_x_continuous(breaks = integer_breaks(10)) +
  scale_y_log10()


library(cowplot)

# Statistical results for flakes. (a-d) The counts of flakes for different mass, different length/thickness ratios (dotted line shows the median value), cortex proportion and number of dorsal scars. (e-g) Density distribution of flakes for different maximum dimension, thickness, and width. (h) Comparison of density distributions of flakes with and without faceted platforms. (i) Box plots showing the mass difference between flakes by scar number. 

figure_9 <- 
plot_grid(
  flake_hist_mass,    flake_lw_ratio, flake_cortex_prop, 
  flake_scar_count,   flake_dens_max, flake_thick_dens,
  flake_width_dens,   flake_dens_max_fac,  flake_scar_by_mass,
  nrow = 3,
  ncol = 3,
  labels = "auto"
)

save_plot(here::here("analysis/figures/figure_9.png"), 
          figure_9,
          base_aspect_ratio = 1,
          base_height = 9,
          base_width = 9)

knitr::include_graphics(here::here("analysis/figures/figure_9.png")) 

Comparing Levallois, Quina, and Non-Levallois flakes

# Note the path that we need to use to access our data files when rendering this document
library(readxl)
library(tidyverse)

file_name <-
  here::here("analysis/data/raw_data/artefacts-after Sam_Dec2018_quina_modified-2019.xls")

flakes <- as.data.frame(read_excel(file_name, sheet = "flake basics"))
cores <- as.data.frame(read_excel(file_name, sheet = "core basics"))
debris <- as.data.frame(read_excel(file_name, sheet = "chunk&debris"))
retouch <- as.data.frame(read_excel(file_name, sheet = "retouch"))
core_flake_scars_1 <- as.data.frame(read_excel(file_name, sheet = "core flake scars"))
edge_angle_1 <- as.data.frame(read_excel(file_name, sheet = "edge angle"))
retouched_pebble_chunk <- as.data.frame(read_excel(file_name, sheet = "Pebble&chunk"))

scar_dir <- as.data.frame(read_excel(file_name, sheet = 'scar direction'))

core_flake_scars <- core_flake_scars_1[core_flake_scars_1$number %in% cores$number, ]

edge_angle <- edge_angle_1[edge_angle_1$number %in% retouch$number, ]

Prepare and clean the data

# prepare the data --------------------------------------------------

require(Ckmeans.1d.dp)
result <- Ckmeans.1d.dp(flakes$mass[!is.na(flakes$mass)])

n_clusters <- length(result$size)
n_clusters_30 <- length(result$size[result$size > 30])

# assign cluster ID to table
flakes_clustered <-
  flakes %>%
  filter(!is.na(mass)) %>%
  mutate(cluster = result$cluster) %>%
  right_join(flakes) %>%
  filter(cluster %in% 1:n_clusters_30)

# clean data to help with analysis
flakes$retouched <- ifelse(grepl("ret", flakes$type1), "retouched", "unretouched")
flakes <- flakes[!flakes$material %in% c('sandstone', 'quartz'),]

Separate Levallois and Non-Levallois flakes for comparison

library(tidyverse)
library(ggbeeswarm)
library(papaja)

# clean data to help with analysis
flakes$retouched <- ifelse(grepl("ret", flakes$type1), "retouched", "unretouched")
flakes <-
  flakes[!flakes$material %in% c('sandstone', 'quartz'), ] %>%
  mutate(complet_flk = ifelse(
    type1 %in% c(
      "flake",
      "blade",
      "KBW",
      "tablet",
      "leva",
      "leva?",
      "DBD",
      "crest"
    ),
    "complete flake",
    ifelse(
      type1 %in% c("flake brk", "blade brk", "LVFB", "proximal"),
      "flake breaks",
      "other"
    )
  )) 

# how many levallois?
# unique(flakes$type1)

leva <- flakes %>% 
  filter(grepl("LVT|LVF|LVFB", type1)) %>%
  filter(!grepl("LVF_", type1))

non_leva <- flakes %>% 
  filter(!grepl("LVT|LVF|LVFB", type1))

flakes_L <-  
  mutate(flakes, L = ifelse(grepl("LVT|LVF|LVFB", type1), 
                            "Levallois", "Non-Levallois"))  

Results of descriptive statistics for Levallois and non-Levallois flakes. ‘PLF’ stands for preferential Levallois flake; ‘CF’ stands for complete flake. This is table 1 in the main text

flakes_complete_non_Lev <- 
  flakes_L %>% # only complete flakes for non-Levallois
  filter(complet_flk == "complete flake")

flakes_Lev  <- 
  flakes_L %>% 
  filter(L == "Levallois")

# combine both here
flakes_Lev_non_Lev <- 
bind_rows(flakes_complete_non_Lev,
          flakes_Lev) %>% 
  as_tibble()
  

# check the CVs
CV <- function(the_vector){
  mean_ <- mean(the_vector, na.rm = TRUE)
  sd_ <- sd(the_vector, na.rm = TRUE)
  cv <- (sd_/mean_) * 100
  round(cv, 3)
}

# prepare mean and sd functions
mean_t <- function(x) mean(x, na.rm = TRUE)
sd_t <- function(x) sd(x, na.rm = TRUE)

# for the table 'Results of descriptive statistics for Levallois and non-Levallois flakes. ‘PLF’ stands for preferential Levallois flake; ‘CF’ stands for complete flake.'
compare_lev_non_lev_vars_tbl <- 
flakes_Lev_non_Lev %>% 
  as_tibble() %>% 
  dplyr::select(L, where(is.numeric)) %>% 
  group_by(L) %>% 
  summarise(across(where(is.numeric),
            list(mean = mean_t,
                 sd = sd_t,
                 cv = CV),
            .names = "{.fn}.{.col}")) %>%
  pivot_longer(-L) %>%
  separate(col = name,
           into = c("stat", "var"),
           sep = "\\.") %>% 
  mutate(L = ifelse(L == "Levallois",
                    "PLF",
                    "CF")) %>% 
  unite("lev_stat", c(stat, L), 
        sep = "_") %>% 
  pivot_wider(names_from = lev_stat ) %>% 
  mutate(difference_in_cv = cv_PLF - `cv_CF`) %>% 
  mutate(across(where(is.numeric), round, 2))

library(gt) 
tab <- 
  compare_lev_non_lev_vars_tbl %>%
  gt(rowname_col = "var") %>% 
  tab_spanner(
    label = "Mean (mm)",
    columns = vars(mean_PLF, 
                   mean_CF)) %>% 
  tab_spanner(
    label = "SD",
    columns = vars(sd_PLF, 
                   sd_CF)  
  ) %>% 
  tab_spanner(
    label = "CV (%)",
    columns = vars(cv_PLF, 
                   cv_CF,
                   difference_in_cv) 
  )

tab
Mean (mm) SD CV (%)
mean_PLF mean_CF sd_PLF sd_CF cv_PLF cv_CF difference_in_cv
length 47.43 49.21 18.42 19.97 38.83 40.57 -1.75
max_dimension 55.96 62.80 20.25 24.58 36.18 39.14 -2.96
oriented_width 46.90 50.33 17.91 21.05 38.18 41.82 -3.64
Width_25%max 32.40 36.15 12.40 14.62 38.27 40.44 -2.18
Width_50%max 37.13 41.75 13.37 16.76 35.99 40.15 -4.16
Width_75%max 30.60 34.94 12.49 15.94 40.82 45.63 -4.81
oriented_thickness 12.16 17.96 4.45 9.06 36.64 50.43 -13.79
Thickness_25%max 11.05 16.15 4.44 8.17 40.20 50.58 -10.38
Thickness_50%max 11.95 17.13 4.78 8.71 40.01 50.86 -10.84
Thickness_75%max 10.04 14.31 4.29 7.70 42.74 53.82 -11.09
mass 39.29 70.91 51.86 92.88 132.01 130.98 1.03
platform_width 31.69 33.50 16.38 18.35 51.68 54.77 -3.08
platform_thickness 10.98 13.31 4.66 8.09 42.46 60.83 -18.38
scar_number 3.39 3.23 1.70 1.83 50.23 56.59 -6.35
cortex_percentage 1.93 8.38 5.73 14.75 296.84 176.01 120.84

Are CV values lower for Levallois flakes than for non-Levallois flakes? This result is presented in line in the text in the main text.

compare_lev_non_lev_vars_tbl  %>% 
  dplyr::select(contains("cv_")) %>% 
  stack() %>% 
  wilcox.test(values ~ ind, 
              data = .) %>% 
  apa_print() %>% 
  .$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")
#> [1] "W = 61.00, p = .033"

Are the levallois pieces bigger by mass?

# plotting setting
point_alpha <- 0.1

# is it significant?
t_test_str_l_mass <- 
apa_print(t.test(data = flakes_L, mass ~ L))$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")

# plot
compare_lev_non_lev_mass_plot <- 
ggplot(flakes_L, aes(x = L, 
                     y = mass)) +
  geom_boxplot() +
  geom_quasirandom(alpha = point_alpha) +
  scale_y_log10() + 
  theme_bw() + 
  xlab("") +
  ylab("Mass (g)") +
  ggtitle(paste0("Flake Mass: ", t_test_str_l_mass))

Are the levallois pieces bigger by max dimension?

# is it significant?
t_test_str_l_maxdim <- 
apa_print(t.test(data = flakes_L, 
                 max_dimension ~ L))$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")

# plot
compare_lev_non_lev_maxdim_plot <- 
ggplot(flakes_L, aes(x = L, 
                     y = max_dimension)) +
  geom_boxplot() +
  geom_quasirandom(alpha = point_alpha) +
  scale_y_log10() + 
  theme_bw() + 
  xlab("") +
  ylab("Maximum Dimension (mm)") +
  ggtitle(paste0("Flake Max. Dim.: ", t_test_str_l_maxdim))

Are the levallois pieces bigger by thickness?

# is it significant?
t_test_str_l_thick <- 
apa_print(t.test(data = flakes_L,
                 `Thickness_50%max` ~ L))$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")

# plot
compare_lev_non_lev_thick_plot <- 
ggplot(flakes_L, aes(x = L, 
                     y = `Thickness_50%max`)) +
  geom_boxplot() +
  geom_quasirandom(alpha = point_alpha) +
  scale_y_log10() + 
  theme_bw() + 
  xlab("") +
  ylab("Thickness (mm)") +
  ggtitle(paste0("Flake Thickness: ", t_test_str_l_thick))

Do the Quina tools have different edge angles?

library(reshape2)
edge_angle_2 <- edge_angle

edge_angle_2[, 2:(ncol(edge_angle_2)-1)] <- atan(edge_angle[, 2:(ncol(edge_angle_2)-1)]/2/3)/pi*180*2

df_angle <- melt(edge_angle_2[, 1:(ncol(edge_angle_2)-1)])

df_angle$variable <- substr(df_angle$variable, 1, 9)

colnames(df_angle)[2:3] <- c("section", "angle")

quina_id <- flakes$number[which(flakes$type3 == "quina")]

df_angle$quina <- ifelse(df_angle$number %in% quina_id, "Quina", "Non-Quina")

# siginificance test
t_test_str_q_angle <- 
df_angle %>% 
  drop_na("angle") %>% 
  t.test(data = ., angle ~ quina) %>% 
  apa_print %>% 
  .$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")

compare_qui_non_qui_angle_plot <-
  ggplot(df_angle,
         aes(y = angle,
             x = quina)) +
  geom_boxplot() +
  geom_quasirandom(alpha = point_alpha) +
  theme_bw() +
  ylab("Edge angle (\u00B0)") +
  xlab("") +
  ggtitle(paste0("Flake Edge Angle: ", t_test_str_q_angle) )

Do the Quina tools have different thickness?

df_thick <- melt(flakes[, c("number", "Thickness_25%max", "Thickness_50%max", "Thickness_75%max")])

df_thick$quina <- ifelse(df_thick$number %in% quina_id, "Quina", "Non-Quina")

df_thick$label <- substring(df_thick$variable, 11, 13)

df_thick_tt_test <- 
df_thick %>% 
  nest_by(label) %>% 
  mutate(t_test = list(
    t.test(data = data, value ~ quina) %>% 
      apa_print %>% 
      .$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")
  ))

compare_qui_non_qui_thick_plot <-
  ggplot(df_thick,
         aes(y = value,
             x = quina)) +
  geom_boxplot() +
  geom_quasirandom(alpha = 0.05) +
  geom_text(data = df_thick_tt_test,
            aes(label = t_test,
                x = 1.5, 
                y = 60)) +
  theme_bw() +
  facet_wrap(~label) +
  ylab("Thickness (mm)") + 
  xlab("") 

Do the Quina tools have different GIUR?

# plot GIUR on Quina  --------------------------------------------------

## calculate GIUR
# looking at the GIUR data... there are some duplicate artefact numbers, which is a nuisance

# make the number ID unique
retouch$number <- make.unique(as.character(retouch$number), sep = "_")

retouch_giur_per_section <-
  retouch %>%
  dplyr::select(grep("^number$|_t|_T", names(retouch))) %>%
  gather(variable, value, -number) %>%
  separate(variable, c("section", "tee"), sep = 10) %>%
  pivot_wider(names_from =  tee, values_from = value) %>%
  mutate(t_T = t / T)  %>%
  filter(t_T <= 1)

retouch_giur_per_artefact <-
  retouch_giur_per_section %>%
  group_by(number) %>%
  summarise(mean_giur = mean(t_T, na.rm = TRUE)) %>% 
  mutate(quina = ifelse(number %in% quina_id, "Quina", "Non-Quina"))

df_giur_clustered_tt_test <- 
retouch_giur_per_artefact %>% 
    t.test(data = ., 
           mean_giur ~ quina) %>% 
      apa_print %>% 
      .$statistic %>% 
  str_remove_all("\\$|\\{|\\}|\\\\|mathit")

# code for raincloud plot
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

compare_qui_non_qui_giur_plot <-
  ggplot(retouch_giur_per_artefact) +
  aes(y = mean_giur,
      x = quina) +
  geom_flat_violin(position = position_nudge(x = .25, 
                                             y = 0), 
                   alpha = .8,
                   adjust = 0.35)  +
  geom_quasirandom( 
             size = 1.5, 
             width = 0.25,
             alpha = point_alpha) +
  geom_boxplot(
        width = .25, 
               outlier.shape = NA, 
               alpha = 0.25
               ) +
  ggtitle(paste0("Flake GIUR: ", df_giur_clustered_tt_test)) +
  coord_flip() +
  theme_bw() +
  ylab("GIUR") +
  xlab("")

Plot combining multiple panels, this is Figure 12 is the main text

## plot figures together.
library(cowplot)
figure_leva_quina <-
  cowplot::plot_grid(plot_grid(compare_lev_non_lev_mass_plot, 
                               compare_lev_non_lev_maxdim_plot, 
                               compare_lev_non_lev_thick_plot ,
                               labels = "auto",
                               label_size = 12,
                               ncol = 3),
                     plot_grid( 
                               compare_qui_non_qui_thick_plot ,
                               labels = c("d"),
                               label_size = 12,
                               ncol = 1,
                               rel_widths = c(1,2)),
                     plot_grid(compare_qui_non_qui_angle_plot,
                               compare_qui_non_qui_giur_plot , 
                               labels= c("e", "f"),
                               ncol = 2,
                               #rel_widths = c(1.5, 1),
                               label_size = 12),
                     nrow = 3,
                     align = 'v',
                     axis = 'lr',
                     vjust = 1.5,
                     hjust = -0.5,
                     scale = 0.9) #

save_plot(here::here("analysis/figures/fig-12-Figures-for-Leva-and-quina.png"),
                     figure_leva_quina,
          base_aspect_ratio = 1,
          base_height = 10, #18,
          base_width =  14) #22)

knitr::include_graphics(here::here("analysis/figures/fig-12-Figures-for-Leva-and-quina.png"))

Retouch

This is figure 14 in the main text


###########################################################################

### Results for retouch - Figure 14

file_name <-
  here::here("analysis/data/raw_data/artefacts-after Sam_Dec2018_quina_modified-2019.xls")

flakes <- as.data.frame(read_excel(file_name, sheet = "flake basics"))
cores <- as.data.frame(read_excel(file_name, sheet = "core basics"))
debris <- as.data.frame(read_excel(file_name, sheet = "chunk&debris"))
retouch <- as.data.frame(read_excel(file_name, sheet = "retouch"))
core_flake_scars_1 <- as.data.frame(read_excel(file_name, sheet = "core flake scars"))
edge_angle_1 <- as.data.frame(read_excel(file_name, sheet = "edge angle"))
retouched_pebble_chunk <- as.data.frame(read_excel(file_name, sheet = "Pebble&chunk"))

scar_dir <- as.data.frame(read_excel(file_name, sheet = 'scar direction'))

core_flake_scars <- core_flake_scars_1[core_flake_scars_1$number %in% cores$number, ]

edge_angle <- edge_angle_1[edge_angle_1$number %in% retouch$number, ]



# calculate the edge angle

edge_angle_2 <- edge_angle

edge_angle_2[, 2:(ncol(edge_angle_2)-1)] <- atan(edge_angle[, 2:(ncol(edge_angle_2)-1)]/2/3)/pi*180*2

df_angle <- melt(edge_angle_2[, 1:(ncol(edge_angle_2)-1)])

df_angle$variable <- substr(df_angle$variable, 1, 9)

colnames(df_angle)[2:3] <- c("section", "angle")

# plot the angles

print("edge angles")
#> [1] "edge angles"

p_edge_angle <- ggplot(df_angle,
                       aes(x = angle,
                           col = section)) +
  geom_density() + theme_bw() +
  xlab("Edge angle (\u00B0)") + ylab("Density") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        legend.text = element_text(size = 15),
        legend.position = c(0.85, 0.72),
        legend.title = element_blank())+ theme(plot.margin=margin(10,10,10,10))

## find quina and compare their angles.

quina_id <- flakes$number[which(flakes$type3 == "quina")]

df_angle$quina <- ifelse(df_angle$number %in% quina_id, "quina", "non-quina")

p_edge_quina <- ggplot(df_angle, aes(x = angle, col = quina)) + geom_density() +
  theme_bw() +
  xlab("Edge angle (\u00B0)") + ylab("Density") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        legend.text = element_text(size = 16))+
  theme(legend.position = c(0.85, 0.85),
        legend.title = element_blank())+ theme(plot.margin=margin(10,10,10,10))

## compare thickness of quina and other

df_thick <- melt(flakes[, c("number", "Thickness_25%max", "Thickness_50%max", "Thickness_75%max")])

df_thick$quina <- ifelse(df_thick$number %in% quina_id, "quina", "non_quina")

df_thick$label <- substring(df_thick$variable, 11, 13)

p_thick_quina <- ggplot(df_thick, aes(x = value, col = quina)) + geom_density() +
  theme_bw() +
  facet_wrap(~label) +
  xlab("Thickness (mm)") + ylab("Density") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        legend.text = element_text(size = 15),
        strip.text = element_text(size = 18),
        legend.position = c(0.9, 0.9),
        legend.title = element_blank()) + theme(plot.margin=margin(10,10,10,10))

# plot the thickness of 50% max

p_thick_quina_2 <- ggplot(df_thick[df_thick$variable =="Thickness_50%max", ], aes(quina, value)) +
  geom_boxplot() + theme_bw() +
  xlab("Thickness at 50% max. dim. (mm)") + ylab("Density") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))+ theme(plot.margin=margin(10,10,10,10))

## count the number of shapes of edges.

edge_shapes <- unlist(strsplit(retouch$`edge shape`, ","))

p_edge_shape <- data.frame(table(edge_shapes)) %>%
  ggplot(aes(x = reorder(edge_shapes, Freq), y = Freq)) +
  geom_bar(stat = "identity") +
  ylab("n") +
  xlab("Edge type") +
  theme_bw() +theme(text = element_text(size = 20)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))+
  theme(plot.margin=margin(10,10,10,10))


## number of layers

p_layers <- ggplot(retouch, aes(`number of layers`)) +
  geom_bar() + theme_bw() +
  xlab("Number of layers") +
  ylab("Count") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))+ theme(plot.margin=margin(10,10,10,10))

## number of edges

p_edge_number <-  ggplot(retouch, aes(`number of edge`)) +
  geom_bar() + theme_bw() +
  xlab("Number of edges") +
  ylab("Count") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25))+
  theme(plot.margin=margin(10,10,10,10))

## number of layers for each number of edge

df_layer_edge <- melt(retouch[, 1:3])

retouch$layers <- ifelse(retouch$`number of layers` < 2, "layer = 1", "layer > 1")

retouch$edges <- ifelse(retouch$`number of edge` > 3, "edge no. > 3",
                        paste("edge no. = ", retouch$`number of edge`))

p_edge_layer <- data.frame(table(retouch[,c("layers", "edges")])) %>%
  mutate(edges =  stringr::str_remove(edges, "edge no. ")) %>%
  mutate(edges =  stringr::str_squish(stringr::str_remove(edges, "="))) %>%
  mutate(edges =  factor(edges, levels = c("0", "1", "2", "3", "> 3"))) %>%
  ggplot(aes(x = edges,
             y = Freq,
             fill = layers)) +
  geom_bar(stat = "identity") +
  ylab("n") +
  xlab("Number of edges") +
  theme_bw() +
  theme(axis.text.x = element_text(size = 20),
                    axis.text.y = element_text(size = 20),
                    axis.title.x = element_text(size = 25),
                    axis.title.y = element_text(size = 25)) +
  theme(plot.margin=margin(10,10,10,10)) +
  theme(strip.text = element_text(size = 18),
        legend.position = c(0.8, 0.9),
        legend.title = element_blank())

# scar number and direction

## count how many scars from each direction

## How many % flakes that have dorsal scars from at least 3 complete different directions

# print("How many % flakes that have dorsal scars from at least 3 complete different directions")

# table(unlist(scar_dir[,2:8]))


## calculate GIUR

# text says ". Most specimens were extensively retouched, i.e.,
# more than XX% have a GIUR value greater than 0.X. "

# we need to compute the GIUR as T/t for each sector, then compute the mean of all sectors
# for each artefact

# looking at the GIUR data... there are some duplicate artefact numbers, which is a nuisance

# make the number ID unique
retouch$number <- make.unique(as.character(retouch$number), sep = "_")

retouch_giur_per_section <-
  retouch %>%
  dplyr::select(grep("^number$|_t|_T", names(retouch))) %>%
  gather(variable, value, -number) %>%
  separate(variable, c("section", "tee"), sep = 10) %>%
  pivot_wider(names_from =  tee, values_from = value) %>%
  mutate(t_T = t / T)  %>%
  filter(t_T <= 1)

retouch_giur_per_artefact <-
  retouch_giur_per_section %>%
  group_by(number) %>%
  summarise(mean_giur = mean(t_T, na.rm = TRUE))

# mean(retouch_giur_per_artefact$mean_giur)

# how many higher than some value
# sum(retouch_giur_per_artefact$mean_giur >= 0.5) / nrow(retouch_giur_per_artefact) * 100

# join metrics variables with GIUR data
library(stringr)
retouch_giur_per_artefact_with_metrics <-
  left_join(retouch_giur_per_artefact,
            flakes,
            by = "number") %>%
  # filter out those with repeated IDs
  filter(!str_detect(number, "_"))

# the average max dimension of retouched piece
average_max_dimension_of_retouch <- round(mean(retouch_giur_per_artefact_with_metrics$max.dimension, na.rm = TRUE),2)
# what percentage of retouched pieces have a GIUR greater than N
N <-  0.5
percentage_with_giur_greater_than_N <-
  round(mean(retouch_giur_per_artefact$mean_giur >= N)  * 100, 0)

#----------
# df_giur <- retouch[, c(1,8,10,12,14,16,18,20,22)]

# df_giur[, -1] <- df_giur[, -1] / retouch[, c(9,11,13,15,17,19,21,23)]

df_giur <- retouch_giur_per_artefact_with_metrics

df_giur_clustered <- left_join(df_giur, flakes_clustered[, c("number", "cluster")])

df_giur_clustered_melt <- melt(df_giur_clustered[,-1], id.vars = "cluster")

## compare GIUR for different mass categories

# print("compare GIUR for different mass groups")

df_giur_clustered_melt$label <-
  paste("Cluster",df_giur_clustered_melt$cluster)

df_giur_clustered_melt <-
df_giur_clustered_melt %>%
  mutate(value = parse_number(value)) %>%
  filter(variable == "mean_giur") %>%
  filter(!is.na(cluster))

## calculate the median GIUR for each group

#print("the median GIUR for each group")

#tapply(df_giur_clustered_melt$value, df_giur_clustered_melt$cluster, median, na.rm=TRUE)

df_giur_clustered_melt_medians <-
df_giur_clustered_melt %>%
  group_by(label) %>%
  summarise(median_giur = median(value))

giur_overall_median <- tibble(giur_overall_median = median(df_giur_clustered_melt$value))

library(ggbeeswarm)
p_GIUR_mass <-
  ggplot(na.omit(df_giur_clustered_melt[df_giur_clustered_melt$value <= 1,]),
                      aes(label, value)) +
  geom_boxplot() +
  geom_quasirandom(alpha = 0.2) +
  theme_bw() +
  ylab("GIUR") +
  xlab("") +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25)) + theme(plot.margin=margin(10,10,10,10))

# Invasiveness index
# print("results for invasiveness index")

df_invasive <- retouch[, c(1, 24:39)]

# need to compute the II as sum of each section divided by the number of sections
df_invasive_ii <-
df_invasive %>%
  rowwise() %>%
  mutate(II = sum(section_1_II,
                  section_2_II,
                  section_3_II,
                  is.numeric(section_4_II),
                  section_5_II,
                  section_6_II,
                  section_7_II,
                  section_8_II,
                  section_9_II,
                  section_10_II,
                  section_11_II,
                  section_12_II,
                  section_13_II,
                  section_14_II,
                  section_15_II,
                  section_16_II,
                  na.rm = TRUE) / 16) %>%
  dplyr::select(number, II)

df_invasive_clustered <- left_join(df_invasive_ii, flakes_clustered[, c("number", "cluster")])

df_invasive_clustered_melt <-
  melt(df_invasive_clustered[,-1], id.vars = "cluster") %>%
  filter(!is.na(cluster)) %>%
  mutate(label = paste0("Cluster ", cluster))

II_overall_median <- tibble(II_overall_median = median(df_invasive_clustered_melt$value))

p_II <-
  ggplot(df_invasive_clustered_melt) +
         aes(as.factor(label),
             value) +
  geom_boxplot() +
  geom_quasirandom(alpha = 0.1) +
  theme_bw() +
  xlab("") +
  ylab("II") +
  ylim(0, 0.5) +
  theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25)) + theme(plot.margin=margin(10,10,10,10))


### plot together

figure_retouch <-
  cowplot::plot_grid(plot_grid(
                              #p_flake_ret_dim,
                               p_edge_shape,
                               p_edge_angle, ncol = 3,
                               labels = "auto",
                               label_size = 26,
                               align = 'vh',
                               axis = 'lr'),
                     plot_grid(p_GIUR_mass,
                               p_II,
                               ncol = 1, labels = c("d", "e"), label_size = 26),
                     plot_grid(
                               p_edge_number,
                               p_edge_layer,
                               ncol = 2, rel_widths = c(1.5,1.5),
                               labels = c( "f", "g"),
                               label_size = 26,
                               align = 'vh',
                               axis = 'lr'),
                     ncol = 1,
                     align = 'vh',
                     axis = 'lr',
                     #rel_heights = c(1, 1.5),
                     vjust = 1,
                     hjust = -0.5,
                     scale = 0.9) #

save_plot(here::here("analysis/figures/Figure-14-retouch.png"), figure_retouch,
          base_aspect_ratio = 1,
          base_height = 21,
          base_width = 22)

knitr::include_graphics(here::here("analysis/figures/Figure-14-retouch.png"))

Colophon

This report was generated on 2024-03-13 15:43:59.729765 using the following computational environment and dependencies:

# which R packages and versions?
if ("devtools" %in% installed.packages()) devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.3.3 (2024-02-29)
#>  os       macOS Sonoma 14.1.2
#>  system   x86_64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Asia/Ho_Chi_Minh
#>  date     2024-03-13
#>  pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package       * version    date (UTC) lib source
#>  beeswarm        0.4.0      2021-06-01 [1] CRAN (R 4.3.0)
#>  bslib           0.6.1      2023-11-28 [1] CRAN (R 4.3.0)
#>  cachem          1.0.8      2023-05-01 [1] CRAN (R 4.3.0)
#>  cellranger      1.1.0      2016-07-27 [1] CRAN (R 4.3.0)
#>  Ckmeans.1d.dp * 4.3.5      2023-08-19 [1] CRAN (R 4.3.0)
#>  class           7.3-22     2023-05-03 [1] CRAN (R 4.3.3)
#>  classInt      * 0.4-10     2023-09-05 [1] CRAN (R 4.3.0)
#>  cli             3.6.2      2023-12-11 [1] CRAN (R 4.3.0)
#>  colorspace      2.1-0      2023-01-23 [1] CRAN (R 4.3.0)
#>  cowplot       * 1.1.3      2024-01-22 [1] CRAN (R 4.3.2)
#>  devtools        2.4.5      2022-10-11 [1] CRAN (R 4.3.0)
#>  digest          0.6.34     2024-01-11 [1] CRAN (R 4.3.0)
#>  dplyr         * 1.1.4      2023-11-17 [1] CRAN (R 4.3.0)
#>  e1071           1.7-14     2023-12-06 [1] CRAN (R 4.3.0)
#>  ellipsis        0.3.2      2021-04-29 [1] CRAN (R 4.3.0)
#>  evaluate        0.23       2023-11-01 [1] CRAN (R 4.3.0)
#>  fansi           1.0.6      2023-12-08 [1] CRAN (R 4.3.0)
#>  farver          2.1.1      2022-07-06 [1] CRAN (R 4.3.0)
#>  fastmap         1.1.1      2023-02-24 [1] CRAN (R 4.3.0)
#>  fitdistrplus  * 1.1-11     2023-04-25 [1] CRAN (R 4.3.0)
#>  forcats       * 1.0.0      2023-01-29 [1] CRAN (R 4.3.0)
#>  fs              1.6.3      2023-07-20 [1] CRAN (R 4.3.0)
#>  generics        0.1.3      2022-07-05 [1] CRAN (R 4.3.0)
#>  ggbeeswarm    * 0.7.2      2023-04-29 [1] CRAN (R 4.3.0)
#>  ggplot2       * 3.5.0      2024-02-23 [1] CRAN (R 4.3.2)
#>  glue            1.7.0      2024-01-09 [1] CRAN (R 4.3.0)
#>  gridExtra     * 2.3        2017-09-09 [1] CRAN (R 4.3.0)
#>  gt            * 0.10.1     2024-01-17 [1] CRAN (R 4.3.0)
#>  gtable          0.3.4      2023-08-21 [1] CRAN (R 4.3.0)
#>  here            1.0.1      2020-12-13 [1] CRAN (R 4.3.0)
#>  highr           0.10       2022-12-22 [1] CRAN (R 4.3.0)
#>  hms             1.1.3      2023-03-21 [1] CRAN (R 4.3.0)
#>  htmltools       0.5.7      2023-11-03 [1] CRAN (R 4.3.0)
#>  htmlwidgets     1.6.4      2023-12-06 [1] CRAN (R 4.3.0)
#>  httpuv          1.6.14     2024-01-26 [1] CRAN (R 4.3.2)
#>  jquerylib       0.1.4      2021-04-26 [1] CRAN (R 4.3.0)
#>  jsonlite        1.8.8      2023-12-04 [1] CRAN (R 4.3.0)
#>  KernSmooth      2.23-22    2023-07-10 [1] CRAN (R 4.3.3)
#>  knitr         * 1.45       2023-10-30 [1] CRAN (R 4.3.0)
#>  labeling        0.4.3      2023-08-29 [1] CRAN (R 4.3.0)
#>  later           1.3.2      2023-12-06 [1] CRAN (R 4.3.0)
#>  lattice         0.22-5     2023-10-24 [1] CRAN (R 4.3.3)
#>  lifecycle       1.0.4      2023-11-07 [1] CRAN (R 4.3.0)
#>  lubridate     * 1.9.3      2023-09-27 [1] CRAN (R 4.3.0)
#>  magrittr        2.0.3      2022-03-30 [1] CRAN (R 4.3.0)
#>  MASS          * 7.3-60.0.1 2024-01-13 [1] CRAN (R 4.3.3)
#>  Matrix          1.6-5      2024-01-11 [1] CRAN (R 4.3.3)
#>  memoise         2.0.1      2021-11-26 [1] CRAN (R 4.3.0)
#>  mime            0.12       2021-09-28 [1] CRAN (R 4.3.0)
#>  miniUI          0.1.1.1    2018-05-18 [1] CRAN (R 4.3.0)
#>  munsell         0.5.0      2018-06-12 [1] CRAN (R 4.3.0)
#>  papaja        * 0.1.2      2023-09-29 [1] CRAN (R 4.3.0)
#>  pillar          1.9.0      2023-03-22 [1] CRAN (R 4.3.0)
#>  pkgbuild        1.4.3      2023-12-10 [1] CRAN (R 4.3.0)
#>  pkgconfig       2.0.3      2019-09-22 [1] CRAN (R 4.3.0)
#>  pkgload         1.3.4      2024-01-16 [1] CRAN (R 4.3.0)
#>  plyr            1.8.9      2023-10-02 [1] CRAN (R 4.3.0)
#>  png             0.1-8      2022-11-29 [1] CRAN (R 4.3.0)
#>  profvis         0.3.8      2023-05-02 [1] CRAN (R 4.3.0)
#>  promises        1.2.1      2023-08-10 [1] CRAN (R 4.3.0)
#>  proxy           0.4-27     2022-06-09 [1] CRAN (R 4.3.0)
#>  purrr         * 1.0.2      2023-08-10 [1] CRAN (R 4.3.0)
#>  R6              2.5.1      2021-08-19 [1] CRAN (R 4.3.0)
#>  ragg            1.2.7      2023-12-11 [1] CRAN (R 4.3.0)
#>  rbibutils       2.2.16     2023-10-25 [1] CRAN (R 4.3.0)
#>  Rcpp            1.0.12     2024-01-09 [1] CRAN (R 4.3.0)
#>  Rdpack          2.6        2023-11-08 [1] CRAN (R 4.3.0)
#>  readr         * 2.1.5      2024-01-10 [1] CRAN (R 4.3.0)
#>  readxl        * 1.4.3      2023-07-06 [1] CRAN (R 4.3.0)
#>  remotes         2.4.2.1    2023-07-18 [1] CRAN (R 4.3.0)
#>  reshape2      * 1.4.4      2020-04-09 [1] CRAN (R 4.3.0)
#>  rlang           1.1.3      2024-01-10 [1] CRAN (R 4.3.0)
#>  rmarkdown       2.25       2023-09-18 [1] CRAN (R 4.3.0)
#>  rprojroot       2.0.4      2023-11-05 [1] CRAN (R 4.3.0)
#>  rstudioapi      0.15.0     2023-07-07 [1] CRAN (R 4.3.0)
#>  sass            0.4.8      2023-12-06 [1] CRAN (R 4.3.0)
#>  scales          1.3.0      2023-11-28 [1] CRAN (R 4.3.0)
#>  sessioninfo     1.2.2      2021-12-06 [1] CRAN (R 4.3.0)
#>  shiny           1.8.0      2023-11-17 [1] CRAN (R 4.3.0)
#>  stringi         1.8.3      2023-12-11 [1] CRAN (R 4.3.0)
#>  stringr       * 1.5.1      2023-11-14 [1] CRAN (R 4.3.0)
#>  survival      * 3.5-8      2024-02-14 [1] CRAN (R 4.3.3)
#>  systemfonts     1.0.5      2023-10-09 [1] CRAN (R 4.3.0)
#>  textshaping     0.3.7      2023-10-09 [1] CRAN (R 4.3.0)
#>  tibble        * 3.2.1      2023-03-20 [1] CRAN (R 4.3.0)
#>  tidyr         * 1.3.1      2024-01-24 [1] CRAN (R 4.3.2)
#>  tidyselect      1.2.0      2022-10-10 [1] CRAN (R 4.3.0)
#>  tidyverse     * 2.0.0      2023-02-22 [1] CRAN (R 4.3.0)
#>  timechange      0.3.0      2024-01-18 [1] CRAN (R 4.3.0)
#>  tinylabels    * 0.2.4      2023-09-02 [1] CRAN (R 4.3.0)
#>  tzdb            0.4.0      2023-05-12 [1] CRAN (R 4.3.0)
#>  urlchecker      1.0.1      2021-11-30 [1] CRAN (R 4.3.0)
#>  usethis         2.2.3      2024-02-19 [1] CRAN (R 4.3.2)
#>  utf8            1.2.4      2023-10-22 [1] CRAN (R 4.3.0)
#>  vctrs           0.6.5      2023-12-01 [1] CRAN (R 4.3.0)
#>  vipor           0.4.7      2023-12-18 [1] CRAN (R 4.3.0)
#>  withr           3.0.0      2024-01-16 [1] CRAN (R 4.3.0)
#>  xfun            0.42       2024-02-08 [1] CRAN (R 4.3.2)
#>  xml2            1.3.6      2023-12-04 [1] CRAN (R 4.3.0)
#>  xtable          1.8-4      2019-04-21 [1] CRAN (R 4.3.0)
#>  yaml            2.3.8      2023-12-11 [1] CRAN (R 4.3.0)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

The current Git commit details are:

# what commit is this file at? 
if ("git2r" %in% installed.packages() & git2r::in_repository(path = ".")) git2r::repository(here::here())  
#> Local:    master /Users/bmarwick/Downloads/lithicprodswchinalatemiddlepleistocene
#> Remote:   master @ origin (https://github.com/benmarwick/lithicprodswchinalatemiddlepleistocene)
#> Head:     [9435e45] 2022-02-09: update after resubmission to NEE